Calculation of a Deuterium Double Shock Hugoniot from Ab initio Simulations 
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We calculate the equation of state of dense deuterium with two ab initio simulations techniques, 
path integral Monte Carlo and density functional theory molecular dynamics, in the density range 
of 0.67 < p < 1.60 gem -3 . We derive the double shock Hugoniot and compare with the recent 
laser-driven double shock wave experiments by Mostovych et al. Q. We find excellent agreement 
between the two types of microscopic simulations but a significant discrepancy with the laser-driven 
shock measurements. 



Single-shock laser experiments on liquid deuterium 
have probed the equation of state (EOS) up to 
340 GPa jjj and have measured a significantly higher 
compressibility than predicted by standard models like 
Sesame and by more recent experiments using mag- 
netic drives [|. The laser-d riven experimental findings 
also disagree with results from first principles simula- 
tions {§, | 1 and leave only some chemical models in 
agreement {7^ ||. The EOS is of fundamental importance 
to our understanding of many physical phenomena such 
as the Jovian planets, inertial confinement fusion, and 
pulsed-power produced plasmas. 

In the recent laser-driven double shock (DS) experi- 
ments by Mostovych et al. ||, final pressures of 100-600 
GPa have been reached. The results also disagree with 
the Sesame EOS but agree well with the linear-mixing 
model M, Both models are based on an approximate 
free energy function constructed from known theoretical 
limits, e.g., Sesame uses Saha and Thomas- Fermi-Dirac 
theories for the electronic properties and hard sphere ref- 
erence systems to incorporate atoms and molecules. The 
double shock experiments appear to provide an indepen- 
dent confirmation of both the laser-driven single shock 
results [Q as well as some of the chemical models. Fur- 
thermore the findings were interpreted as an indication 
that the dissociation of molecules makes deuterium more 
compressible than predicted by microscopic simulations. 

In this paper, we compare the Mostovych et al. |j| 
DS results with ab initio simulations from path integral 
Monte Carlo (PIMC) and density functional molecular 
dynamics (DFT-MD). We give details of the calculation 
of the secondary shock hugoniot curve from the simula- 
tion results. In this calculation, we use either the PIMC 
results or combine PIMC and DFT data. Both meth- 
ods yield very similar results. We discuss whether the 
DS results indicate an increased compressibility for the 
primary hugoniot curve and therefore support the sin- 



gle shock experiments, despite the different temperature 
and densities being reached. Furthermore, we discuss 
whether the process of molecular dissociation is impor- 
tant. 

PIMC |l(], [ll| is a first principles simulation method 
that works directly with electrons and protons. Except 
for the problem of fermi statistics, it is an exact solution 
of the many-body quantum problem for a finite system 
in thermodynamic equilibrium. We have treated fermi 
statistics by restricting the paths in the region of positive 
density matrix by using a variational density matrix |l2] ] 
consisting of a set of Gaussian single particle density ma- 
trices; a procedure that leads to lower free energy than 
the restrictions previously employed with free particles. 
The PIMC method quantitatively describes well the phe- 
nomena of electron correlation, formation of atoms and 
molecules, and temperature effects. The effect of fermi 
statistics is only important well below the fermi temper- 
ature, thus we have confidence in the results for tem- 
peratures greater than 20 000 K {§. The PIMC results 
become increasingly accurate for higher temperatures. 

For lower temperatures, where the effects of elec- 
tronic excitations are less important, the DFT results 
become more reliable (relative to PIMC). We performed 
fixed- volume molecular dynamics simulations, employing 
finite-temperature DFT using the Mermin |l3[] functional 
with the electron and ion-kinetic temperatures set equal 
and the generalized gradient approximation jL4| (GGA). 
Our study used the VASP plane-wave pseudopotential 
code, which was developed at the Technical University 
of Vienna |H| . This code implements the Vanderbilt ul- 
trasoft pseudopotential scheme |]l6[ |l7|, and the Perdew- 
Wang 91 parameterization of GGAllJ]. 

In both shock experiments, the initial state of liquid 
deuterium is at a density of po = 0.171 gem -3 at atmo- 
spheric pressure and close to zero temperature (20K). 
In our calculation, we match these conditions by set- 
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ting Eq = — 15.886eV per atom |L8|, P = and using 
= Po /m where m is the deuteron mass. If a shock wave 
with shock velocity u s and particle velocity u p is driven 
through the sample at state (Eq, no, Po), the new state 
on the other side of the shock front {E\,ri\,P\) follows 
from conservation of mass, momentum and energy [fl9| : 
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For a given EOS, u s and u p can be derived from, = 
(,/n and — £77 with £ = (Pi — Po)/nom and n 
I - no/ ni. 

Before discussing the DS experiments, we will briefly 
review the comparison with results for the princi- 
ple Hugoniot from laser shock wave experiments 
with ab initio simulations. The Hugoniot curves from 
PIMC |i| and DFT-MD §, § are in fairly good agree- 
ment but differ substantially from the laser-driven ex- 
perimental results Q, which exhibit a significantly in- 
creased compressibility (up to 6p ) compared to predic- 
tions based on empirical models such as Sesame . The 
differences between the experimental findings and ab ini- 
tio simulations can be characterized by the extra internal 
energy per particle E and pressure P needed to shift the 
PIMC hugoniot curve to obtain agreement with the ex- 
perimental results (see Fig. 1 in Q). The relationship of 
the two shifts follows from Eq. [| 8E = |(1 — ni/no) Spv 
where pv = P/n. The shift needed is 5E = 3 eV per atom 
to the internal energy or Spv = — 2cV to the pressure. 
Both shifts are far outside the ab initio error bars (order 
of 0.3 eV). 

In the DS experiments, a laser is used to drive a 
shock wave through the deuterium sample. This pri- 
mary shock propagates through the sample and then re- 
flects off an aluminum witness plate. This drives a shock 
wave through the aluminum plate and causes a secondary 
shock wave traveling in the opposite direction through 
the already shocked deuterium sample. We label the fi- 
nal DS state of the deuterium as (E 2 ,n 2 ,P 2 ). 

Shock waves in aluminum have been studied inten- 
sively |^0| and the hugoniot is well known. We use the 
linear fit formula, 



uncertainty of the aluminum hugoniot (Eq. ^) is approx- 
imately 1%, which corresponds to a negligible (2 GPa) 
error in the final shock pressure. 

To study the propagation of the secondary shock front 
in deuterium, one preferably goes into the frame of the 
shocked deuterium moving with u Pl (u S2 — u S2 — u pi 
and u P2 = u P2 — u Pl ). In this frame, Eqs. [j]|3| relate the 
primary shock state (Ei,ni, Pi) to the secondary shock 
state (E 2 ,n 2 ,P 2 ). 

To compare with the DS experiments fl, we need to 
obtain the secondary hugoniot for a given EOS. Assum- 
ing a primary shock pressure Pi, we calculate the first 
shock state (Ei,m, Pi) and consequently u Pl . Then we 
determine the remaining variables as a function of P 2 
and u P2 . They need to satisfy two equations: the alu- 
minum state has to be on the aluminum hugoniot and 
the deuterium states 1 and 2 have to satisfy the hugoniot 
relation Eq. ||. To solve the two equations for the two 
unknowns, we use a Newton procedure beginning with 
an initial guess. 
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5.386 kms" 1 if u p < 7.557 km s" 1 
6.8263 kins" 1 otherwise. 



Using the initial aluminum density of pAi = 2.71 gem 3 
and Eq. 0, we can relate the particle velocity u PA1 to the 
pressure in the shocked aluminum pai ■ Since the double 
shocked deuterium material and the shocked aluminum 
are in direct contact, the pressure as well as the particle 
velocity must be constant across the interface, u p2 = u PAl 
and p 2 = pai (impedance match). We assume that the 
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FIG. 1: Phase diagram of deuterium jllj with primary and 
secondary Hugoniot curves. The primary shock velocities are 
indicated on the arrows connecting the different shock states 
from Tab. |. The solid lines are isobars. 

Calculating the DS hugoniot curve in Fig. [j] requires 
an accurate EOS in the temperature and density range 
of 10000 < T < 100000 K and 0.67 < p < 1.4 g cur 3 . 
This overlaps with the region where both PIMC and DFT 
simulations have been applied. 

Our novel EOS results calculated with PIMC and DFT 
simulations for the densities corresponding to r s = 1.75 
and 1.5 are shown in Figs. |^ and |^ combined with earlier 
results for r s = 2.0 pj. The overall agreement of the 
two ab initio methods is quite good. The DFT results, 
more accurate at low temperatures, smoothly join onto 
the PIMC results at higher temperatures for each den- 
sity. Thus we can consider combining both methods in 
different ways in the DS calculation (Tab. |). For high 
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TABLE I: Secondary shock pressure as a function of primary 
shock velocity comparing measurements with simulation re- 
sults. 1 and 2 label the method used for the primary and 
secondary shock states. * indicate our most reliable results 
and * where the DFT EOS || needed to be extrapolated to 
higher temperatures. Statistical error estimate in parenthe- 
ses. 



20000 T 40000 



60000 



u B1 
(kms" 1 ) 


Experi- 
ment^] 
P 2 (GPa) 


PIMC 1 
PIMC 2 
P 2 (GPa) 


PIMC 1 
DFT 2 
P 2 (GPa) 


DFT § 1 
PIMC 2 

P 2 (GPa) 


DFT (5J 1 
DFT 2 
P 2 (GPa) 


20 




200(10) 


193(6)* 


231(13) 


223(1) 


25 


400(80) 


307(8)* 


307(6) 


327(7) 


328(1) 


30 


590(90) 


429(7)* 




439(4) 




35 




570(7)* 




566(3) t 




40 




730(5)* 




709(3) f 





FIG. 2: Internal energy per atom vs temperature calculated 
with PIMC and DFT simulations. The curves for r s = 1.75 
and 1.5 were offset by 2 and 4 eV for clarity. 
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FIG. 3: Pressure vs temperature calculated with PIMC and 
DFT simulations. 



shock velocities, we used the PIMC EOS for both shock 
states. This yields the most reliable secondary Hugo- 
niot curve except for very low shock velocities such as 
u s = 20kms , we can also combine a primary state 
from PIMC with a secondary state from DFT. Alter- 
natively, we can replace the PIMC data with the DFT 
EOS H for the first shock state. 

All different ways of combining the EOS from micro- 
scopic simulations lead to very similar results compared 
to the DS measurements [|| as shown in Tab. | and 
Fig. f|. Ab initio methods predict a secondary shock 
pressure about 25% lower than measured experimentally. 
This is outside of the experimental error bars reported 
in |2lj] (which have increased compared to |)| ) . Only for 
u Sl = 25kms _1 is agreement found. 

Both ab initio methods have statistical and finite size 
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FIG. 4: Comparison of experimental and several theoretical 
double shock Hugoniot curves showing the secondary shock 
pressure as function of the primary shock velocity. 



uncertainties. The statistical errors, larger for PIMC, 
are given in Tab. [i]. Finite size errors are more difficult 
to estimate due to the computational demand. However, 
very small finite-size errors, corresponding to only a 2% 
uncertainty of Pi were determined from other DFT sim- 
ulations U . In addition, PIMC results JTTJ indicate that 
the finite-size dependence can not explain the difference 
with single shock experiments. We estimate that both 
types of uncertainties are significantly smaller than the 
deviations from the experiments. 

The following analysis is based on the 15 individ- 
ual shock measurements, which have been reported and 
shown in Fig. ||. Using a % 2 computation leads to \ 2 = 39 
for the agreement with PIMC simulations (for both first 
and second shocks) and % 2 = 16 for the Ross model M. 
We note that if the experimental error bars were multi- 
plied by 1.5 then one would obtain a reasonable agree- 
ment with the microscopic simulations. As described 
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for the single shock experiments, one can quantify the 
differences between simulation and experiments by es- 
timating shifts SE and Spv. One finds that a shift of 
SE = (4.0 ± 1.5) eV or Spv = -(3 ± l)eV needs to be 
applied to both shock states in order to bring the PIMC 
hugoniot up in pressure to match the measurements. 

We should clarify one important point. All of the ap- 
proaches discussed in the present work (Sesame, PIMC, 
and DFT) contain dissociation within the formulation; 
the latter two include intermediate states and short-lived 
species in the dense fluid at a sophisticated quantum me- 
chanical level. These states are part of the partition func- 
tion and are therefore included in PIMC. The DFT-MD 
method allows one to study the dynamics of dissociation 
and estimate life-times ||. Therefore, the reason for the 
softening of the principal Hugoniot and the behavior of 
the Mostovych et al. double shock experiment cannot 
arise primarily from the dissociation of molecular hydro- 
gen as suggested in || where the comparison of different 
approximate free energy models had lead to this conclu- 
sion. 

In conclusion, we observed significant differences when 
comparing our ab initio simulation results with the av- 
eraged DS measurements in ^J. The discrepancies are 
significantly larger than the combined error bars from the 
experiment and the simulation. As in the case of the sin- 
gle shock experiments, ab initio methods predict results 
which are relatively close to the Sesame model [|| . When 
comparing with the individual shock data points, we find 
agreement within the error bars only at the lower pres- 
sures. The single shock experiments Q suggested that 
deuterium is significantly more compressible than pre- 
dicted by ab initio simulations. A substantial increase in 
internal energy or decrease in pressure would be neces- 
sary to bring the ab initio EOS in agreement with these 
measurements. We find that comparable shifts in energy 
or pressure are required to reach agreement between the 
DS experiments |9(] and microscopic simulation results. 
Therefore, it can be concluded that both types of laser- 
driven experiments are in relative agreement with each 
other and disagree with ab initio methods. Clearly, more 
theoretical and experimental work is necessary to resolve 
the discussed discrepancy. 

Work performed under the auspices of the U.S. Depart- 



ment of Energy though the University of California at the 
Lawrence Livermore National Laboratory (W-7405-Eng- 
48, CSAR subcontract B341494) and at the Los Alamos 
National Laboratory (W-7405-Eng-36). We acknowledge 
useful discussions with Dr. M. Knudson. 



[1] I. B. Da Silva et al, Phys. Rev. Lett. 78, 783 (1997); 

G. W. Collins et al, Science 281, 1178 (1998). 
[2] G. I. Kerley, in Molecular Based Study of Fluids, edited 

by J. M. Hail and G. A. Mansoori (ACS, Washington 

DC, 1983), p. 107. 
[3] M.D. Knudson et. al., Phys. Rev. Lett, (in press). 
[4] B. Militzer and D. M. Ceperley, Phys. Rev. Lett. 85, 

1890 (2000). 

[5] T. Lenosky et al, Phys. Rev. B 61, 1 (2000); L. Collins 

et al, Phys. Rev. 63, 184110 (2001). 
[6] G. Galli et al, Phys. Rev. B 61, 909 (2000); S. Bagnier 

et al, Phys. Rev. E 6302, 5301 (2001). 
[7] M. Ross, Phys. Rev. B 58, 669 (1998). 
[8] D. Saumon and G. Chabrier, Phys. Rev. A 46, 2084 

(1992). 

[9] A. N. Mostovych et al, Phys. Rev. Lett. 85, 3870 (2000). 
[10] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995); D. M. 

Ceperley, J. Stat. Phys. 63, 1237 (1991). 
[11] B. Militzer and D. M. Ceperley, Phys. Rev. E 63, 066404 
(2001). 

[12] B. Militzer and E. L. Pollock, Phys. Rev. E 61, 3470 
(2000). 

[13] N.D. Mermin, Phys. Rev. 137A (1965) 1441. 

[14] J. P. Perdew, in Electronic Structure of Solids, ed. by F. 

Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991). 
[15] G. Kresse and J. Hafner, Phys. Rev. B 47, RC558 (1993); 

G. Kresse and J. Furthrmiller, Comput. Mat. Sci. 6, 15- 

50 (1996); G. Kresse and J. Furthrmiller, Phys. Rev. B 

54, 11169 (1996). 
[16] D. Vanderbilt, Phys. Rev. B 41 7892 (1990). 
[17] G. Kresse and J. Hafner, J. Phys. Condens. Matter 6, 

8245 (1994). 

[18] W. Kolos and L. Wolniewicz, J. Chem. Phys. 41, 3674 
(1964). 

[19] Y. B. Zeldovich and Y. P. Raizer, Physics of Shock Waves 
and High- Temperature Hydrodynamic Phenomena (Aca- 
demic Press, New York, 1966). 
[20] A. Mitchell and W. Nellis, J. Appl. Phys. 52, 3363 (1981) 
[21] A. Mostovych et al, Phys. Plasmas 8, 2281 (2001). 



